eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_approx.mat')]); 

% Figure 2a
bins = 50;    
eta=.1;
e_f = (RGDP_doubledef_hat_fj-1)/eta*10;
fig1 = strcat('Figures/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Fig2a.eps');
figure(1)
histogram(e_f,bins,'Normalization','Probability')
xlabel('$d \ln Y_{f,n}^F$','Interpreter','latex')
ylabel('Frequency','Interpreter','latex')
saveas(gcf,fig1,'epsc2')

% Figure 2 b
clear
rhoT = '3';
etaT = '1.000001';
lambdaT = '1.000001';
alphaT = 4;
varphi = '3';
sigmaT = '1.5';
shockT='WORLDprod_p10';

eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_approx.mat')]); 
eta=.1;
e_Y = (RGDP_doubledef_hat_n(15)-1)/eta*10;
e_f = (RGDP_doubledef_hat_fj-1)/eta*10;
w_f = VA_fj_0/sum(VA_fj_0);
[q,ed] = discretize(log(w_f),20);
for i = 1:20
    ee = e_f(q==i);
    me = mean(ee);
    mbin_e_f(i) = me;
end
fig2 = strcat('Figures/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Fig2b.eps');
figure(2)
plot(ed(2:end),mbin_e_f,'bo','LineWidth',2)
set(gca,'xscale','log')
hline = refline([0 e_Y]);
hline.Color = 'r';
xlabel('$\omega_{f,n}$ bin','Interpreter','latex')
ylabel('Mean $d\ln Y^F_{f,n}$','Interpreter','latex')
saveas(gcf,fig2,'epsc2')